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Abstract 

A technique is presented for the analysis of measured data obtained from a rotating 
microphone rake system. The system is designed to measure the interaction modes of 
ducted fans. A Fourier analysis of the data from the rotating system results in a set of 
circumferential mode levels at each radial location of a microphone inside the duct. Radial 
basis functions are then least-squares fit to this data to obtain the radial mode amplitudes. 
For ducts with soft walls and mean flow, the radial basis functions must be numerically 
computed. The linear companion matrix method is used to obtain both the eigenvalues of 
interest, without an initial guess, and the radial basis functions. The governing equations 
allow for the mean flow to have a boundary layer at the wall. In addition, a nonlinear least- 
squares method is used to adjust the wall impedance to best fit the data in an attempt 
to use the rotating system as an in-duct wall impedance measurement tool. Simulated 
and measured data are used to show the effects of wall impedance and mean flow on the 
computed results. 


I. Introduction 

A diagnostic tool was needed to measure the interaction modes of ducted fans in order to develop fan 
tone noise reduction techniques. Cicon et al. 1 proposed a concept for a continuously rotating microphone 
rake that would be able to obtain all modal amplitudes in a reasonable data acquisition time. The key 
concept is that by slowly rotating the microphone rake in a manner locked to the fan shaft, a Doppler 
shift is imparted to the duct spinning modes that is uniquely based on the mode physics. The resulting 
rotating microphone rake system developed by the NASA Glenn Research Center (called the Rotating Rake 
Mode Measurement System or just the Rotating Rake) has made several advances in the understanding of 
turbofan mode generation. 2 Hall et al. 3 developed the original data analysis scheme for the Rotating Rake 
using closed form solutions that were readily obtained for hard wall, circular ducts with no mean flow. The 
Rotating Rake has been very successful in measuring the modal content propagating in hard wall ducts. The 
data analysis was extended by Sutliff 2 in an attempt to include the effects of a constant mean flow in a soft 
wall duct, whether circular or annular in cross-section. For a constant area duct, the basis functions were 
again obtained using closed form analytical solutions. 

The initial processing of the data obtained from the Rotating Rake follows two basic steps: (1) The 
signals streaming from the microphones in the rake are digitized and a Fourier transform is applied. The 
Fourier coefficients represent the complex amplitudes of the circumferential modes of the measured acoustic 
pressure field in the duct. (2) Each circumferential mode is further decomposed into a series of radial basis 
or mode functions. This series of radial basis functions takes the form of a weighted sum of Bessel functions 
that are the closed formed solution to the governing wave equation in a circular or annular duct for the 
following conditions: (a) the walls of the duct are hard and there is no mean flow, (b) the walls of the duct 
are hard and there is a uniform mean flow, and (c) the walls of the duct are lined with a locally reacting 
impedance and there is no mean flow. For these cases, the orthogonal properties of the Bessel functions 
allow the complex amplitudes of the radial basis functions to be easily obtained from an integration of the 
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pressure disturbance across the area of the duct. These radial basis functions are linearly combined to fit 
the measured circumferential mode data in a least-squares sense. In addition, there is a closed form relation 
between the radial eigenvalue that is part of the argument of the Bessel function and the axial wave number, 
which determines the propagation and/or decay of the mode in the axial direction of the duct. The data 
analysis for constant mean flow in a soft wall duct proposed by Sutliff 2 simplified the analysis by applying the 
orthogonality principle to the radial basis functions in that case, even though the principle did not strictly 
apply. 

For the general duct flow case where a boundary layer exists at the wall whether the wall is hard or soft, 
the previous data analysis using closed form solutions may not be sufficiently accurate. Example radial basis 
functions are shown in Figure 1 for hard wall and soft wall boundary conditions where the flow is varied 
from no flow to plug flow to using a flow profile. The no-flow and plug-flow cases have identical radial basis 
functions for the hard-wall boundary condition and small changes are seen in the radial basis functions for 
the soft wall case. For these soft-wall conditions, it is possible to use the no-flow radial basis functions for 
fitting rotating rake measured data in a soft wall duct within some error, hopefully small. When a mean flow 
profile is included when computing the radial basis functions (See Figure 9 for plots of the flow profile.), the 
results in Figure 1 show a definite change especially for the soft wall conditions. The MA V g = —0.25 profile 
results represent approach conditions for a jet engine of the size typically used on business jets. The plot 
shows about a 1.4 dB peak mode amplitude change from the non-flow-profile case which would then effect the 
level of acoustic power contained in that mode. The desire to be accurate in determining the mode power in 
a soft wall duct with flow gives sufficient reason to further study the effects of shear flow on the data analysis 
where the governing equations do not have closed form solutions. A more general approach is proposed in 
this paper to the data analysis technique for the rotating rake system by developing a numerical approach 
to obtain the radial basis functions. The technique includes wall impedance boundary conditions and allows 
shear flow conditions to exist in the duct. In addition, the analysis provides for the possibility that the wall 
impedance is unknown in the least-square fit of the basis functions to the data. Since the impedance affects 
all the basis functions nonlinearly, the least-squares fit, in this case, is performed iteratively. 



Figure 1. Comparison of radial basis functions computed using m = 4, n = 1, rj = 10.8 with no flow, plug flow, and flow 
profile conditions, (a) Hard wall, (b) Soft wall with Z = 1.7 — 2.2 i. Flow profiles used for calculations are shown in 
Figure 9. 


In the next section, we present the governing equations for the propagation of sound in a cylindrical or 
annular duct while allowing for the presence of a general parallel mean flow. Using both hard and soft wall 
boundary conditions, an eigenvalue problem is formed that when solved produces both the complex axial 
wave numbers, the eigenvalues, and their associated radial basis functions, the eigen- functions. Since the 
eigenvalue appears nonlinearly in the governing equation, we use the linear companion matrix method to 
directly solve the eigenvalue problem. The method is globally convergent with no need for any initial guess 
eigenvalues unlike alternative approaches that use local iteration. Example results show the utility of the 
method. This is followed by a description of the least-squares procedures for solving for both the radial mode 
amplitudes and the wall impedance. The procedures are designed to accommodate the changing behavior 
of mode shapes and radial-mode order as the wall impedance changes. Finally, results are shown for fitting 
numerically generated radial basis functions to both simulated and measured data. For the latter, both plug 
flow and mean profile flow are used in the calculations. 
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II. Numerical Solution 


The goal of analyzing the Rotating Rake data is to obtain the amplitude of the measured modes in the 
duct and, if desired, deduce the impedance of the installed duct wall acoustic treatment. Once the mode 
amplitudes are found, the acoustic power in each mode can be computed and compared to far field radiated 
power levels. For that to occur, two basic processes are followed: 1) determine the radial basis functions by 
numerically solving the eigenvalue problem presented by the governing equation and 2) determine the mode 
amplitudes by a least-squares fit and if the impedance is to be found, use an iteration approach. 


II. A. Eigenvalue Problem 

The numerical analysis of the eigenvalue problem is based on a discretization of the third order equation 
for the acoustic pressure p derived from combining the linearized momentum and energy equations, and 
assuming the static pressure is constant. 


D 3 P 2 L 9{cM) d 2 p 
Dt 3 C dr dx dr 


1 dc 2 D dp 
c 2 dr Dt dr 


D ^2 ' 


= o 


(i) 


where 


D 

Dt 


d 9 

dt + cM dx 


and V 2 is the Laplace operator in cylindrical coordinates. This equation accounts for the compressibility 
of the flow and allows both the Mach number M and the speed of sound c to be functions of the radial 
direction. For subsonic flow in a duct, the conditions are such that the speed of sound is about constant. 
Hence, equation (1) reduces to that given by Eversman 4 (equation 13) in cylindrical coordinates. 

For a circular or annular duct, the soft wall boundary condition for a locally reacting impedance is found 
by using the continuity of particle displacement. The general form of the boundary condition in terms of 
acoustic pressure at r = R 0l the outer wall, is 


d 2 p 
dt dr 


-~pK 


9 9 

dt + ° M dx 


(2) 


where A* w is the outer wall admittance and p is the mean flow density. (The wall impedance is related to 
the the wall admittance by Z* = 1/A*.) This condition reduces to the pure admittance boundary condition 
when M = 0 at the wall for a boundary layer flow and reduces to the hard wall boundary condition when 
A* ow = 0. For a circular duct, the boundary condition at the duct centerline, r = 0, is given by 



p = 0 


to = 0 
to 0 


(3) 


where m is the circumferential mode number. An annular duct may also have a soft wall boundary condition 
at the inner wall, r = Ri. 


d 2 p 

dtdr 


PA* U 


9 9 

di +cM d. x 


(4) 


The sign difference between equations (2) and (4) follows from the convention that in the definition of the 
wall impedance, the particle velocity is defined as positive into the locally reacting impedance. 

The analysis of equation (1) starts by using the separation of variables solution form for the acoustic 
pressure 

p(x, 9, r, t) = PmnEmn(r)e i (“ t+me - k •*) (5) 


where we let E mn be a general function of r that forms the basis function for the least-squares fit, u> is the 
radian frequency, and k x is the axial wave number which also depends on m and n. Equation (5) represents 
one term in a double summation over all circumferential mode numbers to and radial mode numbers n. Since 
equation (1) is linear, it can be solved using one mode at a time. Substitute into equation (1) and divide 
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out the exponential functions and the amplitude coefficient to obtain an equation that is used to determine 
the shape of the radial basis function E for each pair of (m, n) mode numbers. 


1 3 


1 - cM ■ 


E 


k x d(cM) dE 


1 - cM ■ 


dr dr 


d 2 E 


dr 2 


1 dc 2 
c 2 dr 
1 dE 
r dr 


1 - cM ■ 


dE 

dr 


-^E-klE 


= 0 


( 6 ) 


With boundary conditions defined by equations (2) to (4), the resulting equation becomes an eigenvalue 
problem where the eigenvalue is the complex axial wave number and the eigen-functions no longer have 
the form of Bessel functions that were found for the previous closed form solutions. Using finite difference 
approximations, equation (6) is discretized using N + 1 nodes resulting in a matrix equation that, including 
the appropriate boundary conditions, is written as 


D 3 (fc x )E = 0 


( 7 ) 


where is a lambda matrix of degree three representing a scalar polynomial with matrix coefficients 


D 3 (k x ) = C 0 fc 3 + Cr k 2 x + C 2 k x + C 3 


(8) 


and E is the vector of acoustic pressure basis function values at each node point. 

The eigenvalues of equation (7) are easily found using the linear companion matrix method. Following 
Bridges & Morris, we assemble the companion matrix. The C matrices in equation (8) are populated by 
values obtained from the discretization of equation (6) and its boundary conditions. The computational grid 
has uniform grid spacing after any necessary grid stretching in the physical domain to allow more grid points 
within the boundary layer, if it exists, of a flow. Since Co multiplies the highest powered eigenvalue k x , it is 
found that the boundary points give a zero value on the main diagonal of Co resulting in a singular matrix. 
To avoid the singular matrix and allow the technique to succeed, a transformation of the form 

A = (zi + k x )/(z 2 - k x ) 

is applied to equation (8) leading to the modified equation 

[AqA 3 + AiA 2 + A-2 A + A 3 ] E = 0 (9) 


where 


Ao = Co~f T Ci^2 T C 2 Z 2 T C 3 
Ai = — 3CoZi~| + Ci (z 2 — 2 Z 1 Z 2 ) T C2(2z2 — .zi) + 3C 3 
A 3 = 3 Co~i^ 2 + Ci(^i — 2z±Z2) + C2(^2 — 2zi) + 3C 3 
A 3 = -C 0 zf + C lZ \ - C 2 *1 + C 3 


In the transformation, Z 2 cannot be equal to an eigenvalue k x . Otherwise, we use the values z\ = Z 2 = 1. 
With this transformation, the Ao matrix is no longer singular. We can now left multiply equation (9) by 
Aq 1 , define the vectors Ei = AE and E 2 = AEi, and rearrange the system of equations to obtain 


[A - AI] E = 0 


where 


A = 


— Aq 'Ai 


I 

0 


— Ao'A 2 -Ao'A 3 

0 0 

I 0 


is the companion matrix of order 3 N + 3 and 


E = 


E 2 


Ei 



(10) 
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is the eigen-vector for the companion matrix. The axial eigenvalues k x and their associated eigen-functions 
are found from using the QR algorithm on the companion matrix A to obtain A and applying the inverse 
transform. The eigen-functions E are obtained from the last third of the eigen-vector E and renormalized 
such that the L 2 -norm of E equals 1. More specifically, the eigenvalues k x are written as k xmn where m 
is the circumferential mode number that is fixed for each solution of equation (6). The eigenvalues, for a 
fixed m, are ordered by the index n representing the radial mode order of the eigenvalue. Thus for a given 
mean flow profile, frequency, wall admittances, and circumferential mode number, we obtain all the axial 
eigenvalues and eigen-functions for each radial mode in the duct to some order of the discretization. The 
eigenvalues for higher modes become very inaccurate. The axial eigenvalues with the smaller magnitudes 
are of the most interest since they are associated with the lower order radial eigen-functions for which the 
Rotating Rake is designed to measure. 

The eigenvalues produced by the companion matrix method are associated with all the modes that can 
propagate and/or decay in both directions in the duct. Since the current technique assumes that modes can 
propagate and decay only in one direction, the appropriate eigenvalues must be identified from the global 
set produced. Using the sign convention used for the exponential functions in equation (5), the following 
definitions are used to identify the appropriate eigenvalues. If k x is complex, then the eigenvalues of physical 
interest are those that decay in the direction of propagation. These are identified by evaluating the imaginary 
part of k x using the conditions 


SjUc} < 0 decay in +x direction 

Sjfca,} > 0 decay in — x direction (11) 

If k x is real, then the group velocity c g = duj/dk x is computed to indicate the direction of propagation. 

c g > 0 propagation in +x direction 

c g < 0 propagation in — x direction (12) 

Using c g to determine the propagation direction takes into account the effect that a mean flow has on the 
propagation. A subset of lower-ordered eigenvalues representing modes propagating and/or decaying in the 
-hr direction is identified from the global set obtained and stored for use in the least-squares process. 


Mode 

Vo et al. 6 

N = 200 


Mode 

Vo et al. 6 

N = 200 

1 

-1.5459i 

-1.5459i 


1 

1.880 - 0.845i 

1.8797 - 0.8450i 

2 

-5.2368i 

-5.2368i 


2 

0.278 - 3.414i 

0.2783- 3.414U 

3 

-8.4775i 

-8.4775i 


3 

0.145- 6.798i 

0.1451 - 6.7979i 

4 

-11.6632i 

-11.6632i 


4 

0.099 - 10.025i 

0.0991 - 10.0248i 

5 

-14.8299i 

-14.8299i 


5 

0.075- 13.21 li 

0.0754- 13.2106i 

6 

-17.9877i 

-17.9877i 


6 

0.061 - 16.379i 

0.0609 - 16.3793i 

7 

-21.1407i 

-21 . 1407i 


7 

0.051 - 19.539i 

0.0511 - 19.5392i 

8 

-24.2908i 

-24.2907i 


8 

0.044 - 22.694i 

0.0440 - 22.6940i 

9 

-27.4388i 

-27.4387i 


9 

0.039 - 25.846i 

0.0387 - 25.8456i 

10 

-30.5856i 

-30.5853i 


10 

0.034 - 28.995i 

0.0345 - 28.9949i 


Hard wall, m = 1, ry = 1.0, M = 0.0 Soft wall, A ow = 1.0 + 1.0*, m = 0, 77 = 1.0, M = 0.0 

Table 1. Comparison of k x eigenvalues for modes Table 2. Comparison of k x eigenvalues for modes propa- 

propagating in the +x direction. gating in the +x direction. 


II. B. Example Eigenvalue Problem Solutions 

The computation of the eigenvalues and eigen-functions using the linear companion matrix method was 
performed with 5-point discretization at all interior and boundary points. This gives a fourth-order accurate 
solution. To validate the solution, test cases were run with increasing complexity. First, the eigenvalues 
were computed for a hard wall, cylindrical duct with no mean flow. The results are shown in Table (1) for 
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the to = 1 circumferential mode with 77 = 1.0, where 77 = ujR 0 /c r is a nondimensional frequency and c r is 
a reference speed of sound. The comparison results are given from the fourth-order Runge-Kutta solution 
of Vo & Eversman. 6 The eigenvalues are identical to at least 3 decimal places. Next, the outer wall is 
made soft with a specific acoustic admittance A ow = p r c r A* w = 1.0 + 1.0* where p r is a reference density. 
The computed eigenvalues are identical to those from Vo & Eversman as shown in Table (2) for the m = 0 
circumferential mode. 

The problem complexity is increased by adding a constant mean flow to the soft wall circular duct. To 
represent the case of flow entering an inlet, the x-axis of the problem points out of the inlet while the flow 
enters in the opposite direction. Hence, the flow Mach number is set at a negative value, M = —0.5. The 
eigenvalues are shown in Table (3). The results are again identical to those from Vo & Eversman. (Note 
that the 10th mode from Vo & Eversman is actually the 14th mode. This was verified by observing the 
eigen- functions . ) 


-(-Mode 

Vo et al. 6 

N = 200 

-Mode 

Vo et al. 6 

N = 200 

1 

0.620- 5.014i 

0.6195- 5.0139i 

1 

0.410 

+ 

1.290i 

0.4101 + 1.2904i 

2 

-5.820 - 3.897i 

-5.8195- 3.8968i 

2 

1.259 

+ 

6.085i 

1.2595 + 6.0852i 

3 

0.445 - 9.187i 

0.4453- 9. 1868i 

3 

1.146 

+ 

9.668i 

1.1457 + 9. 6679i 

4 

0.453- 13.062i 

0.4539- 13.0615i 

4 

1.022 

+ 

13.315i 

1.0218 + 13.3150i 

5 

0.480 - 16.822i 

0.4795 - 16.8217i 

5 

0.943 

+ 

16.977i 

0.9425 + 16.9768i 

6 

0.503- 20.531i 

0.5029 - 20.5307i 

6 

0.891 

+ 

20.635i 

0.8907 + 20.6353i 

7 

0.522 - 24.213i 

0.5220- 24.2131i 

7 

0.855 

+ 

to 

to 

00 

00 

0.8549 + 24.2884i 

8 

0.538 - 27.880i 

0.5376 - 27.8799i 

8 

0.829 

+ 

27.937i 

0.8288 + 27.9368i 

9 

0.550- 31. 537i 

0.5503- 31.5370i 

9 

0.809 

+ 

31.581i 

0.8089 + 31.5814i 

10 

0.589 - 49.754i 

0.5608- 35. 1875i 

10 

0.755 

+ 

49.772i 

0.7933 + 35.2232i 

14 


0.5891 - 49.7539i 

14 




0.7546 + 49.7717i 


Soft wall, A ow = 0.72 + 0.42z, m = 2, rj = 1.0, M = —0.5 


Table 3. Comparison of k x eigenvalues. -|-Mode indicates mode propagating in the -\-x direction against the flow. -Mode 
indicates mode propagating in the -x direction with the flow. 


The possibility of including a shear flow in a duct is the advantage of the above numerical solution 
compared to the previous closed-formed, analytical solutions. This case was tested next. Previous eigenvalue 
results were available in Kousen 7 for a hard wall duct with a mean flow having a 1 /7th power law profile. 

M(r) = M{ 1 - r) (1/7) 

The computed eigenvalues are shown in Table (4) compared to the results from Kousen. The real eigenvalues 
and the imaginary parts of the complex eigenvalues are identical to 3 decimal places. The real parts of the 
complex eigenvalues are the same to 2 to 3 decimal places. 

II. C. Least-Squares Fit 

The microphone data from the Rotating Rake is Fourier analyzed to obtain complex values p m (ri,u>) that 
represent the separation of the acoustic pressure by frequency and circumferential mode order at each radial 
location, ry, of the microphone. These values are expected to be functions of the radial mode content in the 
duct following 


N max - 1 

Pm(h^) ~ ^ ^ A ow , A.^ w1 Oj ) , * = 1 , 2, . . . , iV meas (13) 

n — 0 

Equation (13) represents the attempt to further decompose the measured circumferential mode data into an 
approximate series of weighted radial basis functions computed using the linear companion matrix method 
as discussed above. These functions are not only functions of radial position and frequency, but are also 
functions of the outer duct wall admittance A* w and, in the case of an annular duct, the inner duct wall 


NASA/TM— 2007-2 14966 


6 




Mode 

Kousen 7 

N = 200 

1 

0.81547 

0.81527 

2 

0.76952 

0.76948 

3 

0.72763 

0.72758 

4 

0.65346 

0.65339 

5 

0.54051 

0.54043 

6 

0.36975 

0.36966 

7 

0.06477 

0.06481 

8 

-0.28277 + 0.48650i 

-0.28309 + 0.48624i 

9 

-0.28328 + 0.80439i 

-0.28366 + 0.80420i 

10 

-0.28368 + 1.0539i 

-0.28415 + 1.05367i 

11 

-0.28403 + 1.2756i 

-0.28459 + 1.27546i 


Hard wall, m = 0, r] = 20.0, M = 0.3 


Table 4. Comparison of k x /rj eigenvalues for modes propagating in a duct containing a shear flow. 


admittance A* w . To solve equation (13) for the complex radial mode amplitudes P mn (uj), we use the least- 
squares method. For fixed admittances, the solution is to minimize the residual sum of squares 

s = II P BP || | (14) 


where 


P = 


3?{Pm(?’l)} 


${Pm{rN mea „)} 3{p m (ri)} 


9{Pm(rAT mea J} 


1 T 


is a real vector comprised of the real and imaginary parts of the measured complex circumferential mode 
data of length n p = 2 N meas , 



^{AWo^i)} 

-9{.E m o(ri)} 

'ft { /7,n ( /Y rri 

ax-l)( r l)} 

-%{E m ( Nmax -i)(n)} 

B = 

EmO (r N m eas ) } 

-3{£ m0 (rjv meas )} 

3?{£l m o(ri)} 

^{Pm(N max 

-l)( r iVmeaJ} 

ax-l)( r l)} 

— ' S {E m ( N max —1) { r N mea3 )} 

5t{ E m (N max - 1 ) ( r i ) } 


3{F m0 (nv mea J} 

EmO { r N meas ) } 

E m (N max 

-1 )( r N me as)} 

E m (N m a x — 1) i^Nmeas )} - 


is a real matrix of coefficients derived from the radial basis functions of size n p x n c , n c = 2 N max , and 

r iT 


%{Pm 0 } 


Pm(N max -l )} ~ s {Pm(N max -l)} 


is the real solution vector for the complex radial mode amplitudes of length n c . By ensuring that the number 
of measurement points N meas is greater than the number of radial basis functions N max , the problem is 
overdetermined, n p > n c , and a solution to this linear problem is easily computed using the singular value 
decomposition method. This method provides the best approximation to solving the linear least-squares 
problem for over-determined systems and it is robust, handling problems where B is close to being singular. 
(See Lawson & Hanson 8 for details.) 


II. D. Process for Nonlinear Least-Squares Fitting 

In addition to the process of obtaining the complex modal amplitudes given the radial basis functions, it is 
desirable to deduce the wall admittance from the rotating rake measurements. There are cases where the 
wall admittance is not known or not known precisely where measurements could provide those values. At 
other times, a check of the installed wall admittance is required. To determine the wall admittance from 
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rotating rake data requires more extensive calculations than the linear least-squares fit to determine the 
mode amplitudes. The admittance (or its inverse, the impedance) appears in the boundary conditions of the 
governing equations. As such, the admittance affects the computed mode shapes, which are themselves, in 
essence, power series representations. Thus, the effect of varying the wall admittance appears nonlinearly in 
the mode solution. What was a linear least-squares solution to obtaining the mode amplitudes becomes a 
nonlinear, least-squares iterative solution. However, since only the admittances appear as nonlinear param- 
eters, we can use conditional linearity to separate the fitting of the linear (the amplitude coefficients) and 
the nonlinear parameters. This means that the iterative process required to solve a nonlinear least-squares 
problem is only applied to the nonlinear parameters. 

Let f t be the expectation function defined by the right hand side of equation (13). Using conditional 
linearity, we expand /, into a Taylor’s series about the current admittance values 0 (a vector of length n z ) 
given the current amplitude values P 

f(P,0 + 8)=f(P,0) + B(P,0)8 (15) 


where 

~ df 

B ij = i = l,2,...,n p j = 1,2, ...,n z . (16) 

°Pj P,/3 

In these equations, f(P, (3) is determined from the linear least-squares solution outlined above for the current 
0 resulting in the current P, the best-fit amplitude coefficients for (3. Thus, the amplitude coefficients are 
not held constant, but change to the best-fit values as (3 changes. Since equation (16) must be computed 
using numerical derivatives, additional calculations of f are required at incremental changes of each element 
of /3, all with the appropriate linear least-squares solution for f. Once f and B are computed, the system 
of equations is solved for the change 6 in the admittance values using least squares, computing the residual 
error at the new admittance values, and repeating until the error is minimized. While this iteration process is 
straight forward, in nonlinear least-squares problems it can be unstable and not converge to a minimum error 
solution. For instance, the computed S can be too large, sending f outside the local region of convergence. 
Thus, means are required to control the process. 

Lawson & Hanson 8 detail using the method of Marquardt to solve this nonlinear least-squares problem 
and control the size of the steps that 8 can take. The equations are transformed and amended into the form 


GBD 1 

vn 


(D 8 ) 


G(p-f) 

0 


(17) 


where G is a diagonal matrix used for weighting the data. For example, it is common to give a relatively 
larger weight to terms in (p — f) that are considered more accurate. The diagonal matrix D is computed 
from the square root of the diagonal terms in B T G T GB. This matrix rescales the matrix B to prevent 
erratic behavior in the convergence due to the possibility that B may be nearly singular and causing large 
changes in 8. The £ parameter moderates the step size taken by the impedance parameter [3 from the current 
parameters 0 to the new parameters 0 = 0 + 8. 

Singular value decomposition is again used to solve for 8. Decompose the matrix GBD -1 as 

GBD' 1 = U£V t 


where U contains the left singular vectors, V contains the right singular vectors, and S is a diagonal matrix 
containing the singular values cr*. The step size vector 8 is then solved using the following three equations. 


cr = U t G (p - f) 


n 



i 1 , 2 ,..., n z 


As £ is increased, the step size 8 is reduced to prevent divergence of the iteration. Alternatively, £ can be 
reduced to promote faster convergence when possible. Thus, with appropriate choices for £ as the iteration 


NASA/TM— 2007-2 14966 



progresses, systematic adjustments are made in both the nonlinear admittance parameters and the linear 
amplitude coefficients in order to reduce the error of the fit. Convergence to a solution is achieved when 
either 

s — s 

< e, 


the relative change in the residual sum of square errors, or 


Pi 


< e 


for all i, 


the relative change in all the nonlinear parameters, is less than some chosen convergence tolerance e. 

This procedure for nonlinear least-squares fitting is based on the fact that the measured data is to be fit 
by a series of radial basis functions. As the conditions change in the duct, such as the mean flow and/or 
the wall admittance, the radial basis functions must be recomputed. Oftentimes, one can easily compute 
eigenvalues and eigen-functions by starting with a good initial guess eigenvalue and then tracking it through 
small changes in the duct conditions to the final desired condition. This could be easily accomplished by 
using the shooting or matrix method to solve the eigenvalue problem based on a reduced form of equation 
(6). The inverse power method would then be used to obtain the associated eigen- function. One could easily 
follow one or more modes using this approach with much less calculations then are performed here. However, 
such an approach often fails when the wall admittance is near ‘optimum’ where the theoretical maximum 
attenuation occurs for a particular frequency, geometry and input model content. 10 

0 
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4 

5 

6 
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Figure 2. Trajectories of eigenvalues k x for the 2 lowest radial modes as a function of the wall admittance A ow where 
3{A 0 u) } is held constant. Black: 3{A ou) } = 0.5, Red: } — 0.45, Blue: } — 0.22; m = 2, 77 = 3.6, M — 0.0. The 

region where a double eigenvalue exists is marked with an ‘X’ bounded by = 1.53 and ^{Ao-u,} = 0.45 to 0.50. 

Side plots are radial basis functions at selected k x values. Solid lines are magnitude. Dashed lines are phase. 




The near-optimum admittance case is a test of the numerical methodology where conditions are changing 
rapidly. This is illustrated in Figure 2 where trajectories of the eigenvalue k x are shown as a function of 
the wall admittance. In a manner following Zorumski & Mason, 11 the imaginary part of the admittance 
is held constant along the trajectory and every dot represents a 0.015 increment in the real part from 
the previous dot. Starting from $l{A ow } = 0, the basis functions labeled (2,0) have a real k x , representing 
propagating duct modes. As 5 l{A ow \ increases, the k x trajectories leave the real axis with increasing amounts 
of damping. The basis functions labeled (2, 1) begin with an imaginary k x , indicating that the duct modes are 
cut off. Their eigenvalues k x move toward larger real parts as $l{A ow } increases. The condition for optimum 
admittance corresponds to the location of a double eigenvalue shown in the figure. As can be seen in this 
region, the eigenvalues change rapidly with small changes in admittance. Furthermore, the shapes of the 
basis function convert from one radial order to another as the trajectory passes near the double eigenvalue 
point. 12 In this example, the (2,0) basis functions shown in black convert to (2, 1) basis functions as the 
black trajectory passes on the right side, larger ?R.{k x } values, of the double eigenvalue. Conversely, on the 
left side, lower SR/fez} values, the basis functions convert from a (2, 1) to a (2,0) type of basis function. The 
consequence of this is that the eigenvalues and eigen-functions can change abruptly around near-optimum 
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admittance causing unstable conditions to occur in the nonlinear least-squares iteration. In the fitting 
process, it is the form of the radial basis functions that are fit to the data. It is the linear companion matrix 
method that produces these functions. Since we are not specifically concerned with what order the basis 
functions are listed, we can avoid the potential issue of ‘mode-switching’ by computing both the eigenvalues 
and the eigen-functions at each new admittance condition. The best fit of the basis functions to the data 
is then obtained by starting with one best fit basis function and then including additional best fit basis 
functions one at a time until the desired best fit overall is achieved. This process is illustrated in Figure 3 
where data from measurements taken by the Rotating Rake in a hard wall duct are fit by computed radial 
basis functions. In Figure 3a, the first five radial basis functions are shown in the order computed by the 
linear companion matrix method. The third basis function has the lowest residual sum of square errors, 
s = 197, and is shown fit to the data on the right side with mode amplitude labeled P^o- Next, Figure 3b 
shows the remaining four basis functions. We now fit each of these functions to the data along with the part 
(a) function. The addition of the first basis function provides the smallest error, s = 0.326, and this mode 
amplitude is labeled Pjp. The resulting total fit for the sum of the two modes is shown on the right side. 
Finally, part (c) shows the process repeated for the addition of a third basis function. For the fixed hard-wall 
boundary condition, this illustration is a repeated application of linear least-squares fitting of the computed 
basis functions as each additional basis function is included in the fit equation. In the nonlinear least-squares 
fitting, this process is performed at each iteration step as the wall impedance is varied to achieve best fit. 
While the mode amplitudes that gave minimum error were found directly, finding the best fit impedance 
through iteration requires a reasonable initial starting value. 

III. Numerical Results 

We next show two sets of numerical results for fitting computed radial basis functions to data. First, 
simulated data is used to test the fit process illustrated in Figure 3 for soft wall conditions. Second, data 
taken in the NASA Glenn Advanced Noise Control Fan (ANCF) 13 is used for fitting computed radial basis 
functions. The nonlinear least-squares fit process is used to find the best fit impedance for both constant, 
plug flow conditions and for mean flow profile conditions. 

III. A. Simulated Data 

The Eversman Inlet Radiation code 14 was used to generate mode propagation data along a treated section 
in a cylindrical duct. A (2,0) mode with a given amplitude was specified at the input plane of the duct. 
The duct initially had a hard wall followed by the treated section then reverting to a hard wall before 
the exit of the duct. The code computes the behavior of the initial mode as it propagates through the 
duct including scattering, reflection, and absorption. The magnitudes and phases of the pressures were 
extracted from the numerical data at the duct cross-section equal distance from the two ends of the treated 
section. Nine radial locations were selected to simulate microphone measurement points and the analysis 
was conducted for the m = 2 circumferential mode with r/ = 3.6 and M = 0.0. Two cases were computed: 
one non-optimum impedance, Z = 2.16 — 2.88*, and one near-optimum impedance, Z = 0.6 — 0.2*. Both wall 
impedance conditions used in the Eversman duct code calculations are indicated in Figure 2. Obviously, the 
non-optimum impedance is not near the double eigenvalue condition. 

The duct code numerical data represents a simulation of the processed output from a rotating rake 
measurement, that is the complex pressures for one circumferential mode, the left side of equation (13). To 
illustrate the possibility of using the process for nonlinear least-squares iteration to determine a best fit wall 
impedance for the measured data, least-squares errors were computed as a function of the wall impedance. 
Figure 4 shows contour plots for the fit error using the non-optimum wall impedance data generated by the 
duct code. It is implicit in these results that the mode amplitudes, the P mn (w) factors in equation (13), are 
the best fit values for each wall impedance. At each impedance point used for computation on the grid, the 
process shown in Figure 3 is followed: first for one basis function, then two basis functions, and finally, three 
basis functions. All three cases show the same distinct minimum error region. Hence, these results show 
that the best fit wall impedance should be found in a straight forward manner using an iterative process 
where the error at each step is computed by a linear least-squares code using the most recent impedance. In 
this case, the best fit impedance would only change slightly whether 1, 2, or 3 radial basis functions were 
used in the fit. 
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Mode (4,0) P 4Q = -25.9 + 92. 7i 
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(b) 
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Total Ma § 

Mode (4,0) P 4Q = -26.1 + 92.4i 96.1 

Mode (4,1) P = -56.3- 54.1i 78.0 




Data 

Total 

Mode (4,0) P 40 = -25.2 + 91.8i 
Mode (4,1) P 41 = -55.7 - 54.4i 
Mode (4,2) P 42 = -2.17 + 1.49i 


Mag 

95.2 

77.9 

2.64 


(C) 



Figure 3. Process of fitting radial basis functions F/ mn to data from Rotating Rake measurements in a hard wall duct. 
m = 4, rj = 10.8, and M = —0.115. (a) Best fit one function, (b) Best fit second function given first function, (c) Best fit 
third function given first two functions. All basis functions are real. Mode amplitudes are complex. 
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(a) 1 Radial Basis Function 
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(b) 2 Radial Basis Functions 
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Figure 4. Contours of minimum least-squares fit error as a function of the wall impedance Z for 1, 2, or 3 radial basis 
functions (modes). Data was generated using input mode (2,0), wall impedance Z = 2.16 — 2.88i, 77 = 3.6, and M = 0.0. 


Figure 5 shows the fit error contour plots using the duct code numerical data with the near-optimum 
wall impedance. In this case, more than one basis function is required to obtain a best fit wall impedance 
that has better agreement to the impedance used in the duct code. This result also illustrates that the fit 
process in the vicinity of the near-optimum impedance works well without regard to the ‘mode-switching’ 
behavior shown in Figure 2. The error contours are smooth indicating that the iterative process should not 
have problems progressing through to minimum error. Figure 5a may show the effects of ‘mode-switching’ 
as one basis function is fit in the impedance plane surrounding the double eigenvalue location, but this was 
not studied in detail at this time. For two or more basis functions, the effects of ‘mode-switching’ are not 
noticeable since the two basis functions are always present to balance each other to achieve the best fit. 
Thus, the results in Figures 5b and 5c behave more like the results in Figure 4 where adding more basis 
functions does not change the best fit impedance. 



Re{Z} 

(a) 1 Radial Basis Function (b) 2 Radial Basis Functions (c) 3 Radial Basis Functions 

Figure 5. Contours of minimum least-squares fit error as a function of the wall impedance Z for 1, 2, or 3 radial basis 
functions (modes). Data was generated using input mode (2,0), wall impedance Z — 0.6 — 0.2 i, r) — 3.6, and M — 0.0. 


III.B. Measured Data 

The set of contour plots, shown in Figure 6, are computed fit errors resulting from using measured Rotating 
Rake data taken in the NASA Glenn Advanced Noise Control Fan (ANCF). 13 The wall treatment was a 
linear, single-degree-of-freedom liner with a screen mesh on a 34 percent open area perforate. For a fan speed 
of 1800 RPM, the duct plug flow Mach number was —0.115 and the design wall impedance was Z = 1.7 — 2. 4i 
for a normalized frequency of 77 = 10.8 corresponding to the second harmonic of the blade passage frequency. 
The to = —8 circumferential mode was used in this analysis. Unlike the previous results using computed 
simulated data, these results show significant changes in the error contour plots as the number of radial basis 
functions changes in the fit equation. For one basis function, the minimum error is in the vicinity of the 
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(a) 1 Radial Basis Function (b) 2 Radial Basis Functions (c) 3 Radial Basis Functions 

Figure 6 . Contours of minimum least-squares fit error as a function of the wall impedance Z for 1, 2, or 3 radial basis 
functions (modes). Data from rotating rake measurements, m = — 8 , nominal wall impedance Z = 1.7 — 2.4i, 77 = 10.8, 
and M = —0.115. 


nominal design impedance. Adding a second basis function to the fit gives rise to a second local minimum 
for best fit wall impedance within the search grid. With a third basis function, the two local minimum 
regions are smaller and distinct. While the color contours visually show the regions of minimum error, the 
range of actual error values is decreasing as more basis functions are added to the fit equation. Thus, given 
that each impedance point in Figure 6c has the best fit amplitudes for that impedance, then practically 
any point ‘fits’ the data even though there are two local minimum error regions since all the fit errors are 
reasonably small. To illustrate, the radial basis functions and the best fit modes amplitudes and shapes for 
these two conditions are shown in Figures 7a and 7b. In part (a), the (—8, 0) mode fits the data very well by 
itself. (It should be noted that of the seven microphones the one microphone at r/R = 0.943 was determined 
after the testing to be not working properly. Hence, that data point was weighted such that it was set to 
zero.) Whereas in Figure 7b the addition of the (—8, 1) mode is required to achieve a comparable level of 
fit. This results in a small 8 percent change, 0.7 dB, in the the dominant (—8,0) mode amplitude but the 
best fit impedance has changed about 36 percent. The small changes in mode amplitudes were enough to 
compensate for the changes in the radial basis functions due to the large change in the wall impedance to 
achieve a minimum fit error. 

The same set of measured Rotating Rake data was analyzed again using a measured mean flow profile 
instead of plug flow. The flow profile is shown in Figure 9 with an average Mach number MAvg = —0.115 
(recalling that the flow is in the opposite direction of the positive ic-axis.) The set of fit error contour plots 
are shown in Figure 8 and are contrasted to the equivalent set for plug flow conditions in Figure 6. Clearly 
for this condition the mean flow has affected the conditions for best fit wall impedance. Though it should be 
noted that the range of errors over the grid of impedances again becomes very small when two or more radial 
basis functions are used in the fit. Even though the contour plots can show dramatically different locations 
for best fit impedance, there is actually little change in the basic fit to the data. This is also illustrated in 
Figure 7. By adding only the effect of the mean flow profile at a constant wall impedance, we go from part 
(a) to part (c). The change in the shape of the basis functions is shown and the amplitudes of the modes 
have increased, a 2.5 dB increase for the (—8, 0) mode, but the data is fit quite well. To achieve that last 
10 percent reduction in error to get to the minimum error /best-fit-impedance condition required the real 
part of the impedance to increase 67 percent from 1.42 to 2.37. Thus, while the presence of the mean flow 
profile effects the shape of the radial basis functions, the limited amount of data included in the fit does 
not provide sufficient constraint allowing a large range of amplitude and impedance combinations to provide 
minimal fit errors. Other conditions were run to see if these kinds of changes occurred. The results for best 
fit impedance at three conditions are shown in Figure 10 compared to results from ANCF measured data 
analyzed using closed form basis functions and other measured and computed results previously presented in 
Sutliff. 13 Results are plotted at 747 Hz (at 1400 RPM, M = —0.085, rj = 8.4, and m = —4), 853 Hz (at 1600 
RPM, M = —0.1, 77 = 9.6, and m = —8), and 960 Hz (at 1800 RPM, M = —0.115, r] = 10.8, and m = —8) 
for both the plug flow case and the flow profile case. For the two lower speed cases, the differences in the 
best fit impedance between the plug flow case and the flow profile case is much less than in the highest speed 
case. Also in both cases, the minimum error location was distinct as a function of the wall impedance. It is 
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Figure 7. Changes in radial basis functions and mode fits with changes in wall impedance and flow conditions using 
data from Rotating Rake measurements in a soft wall duct, (a) Plug flow M = —0.115, best fit Z — 1.42 — 2.32 i near 
nominal wall impedance, s = 0.0371. (b) Plug flow M = —0.115, best fit Z — 1.93 — 1.47*, s = 0.0372. (c) Flow profile shown 
in Figure 9 with MAvg = —0.115, wall impedance Z — 1.42 — 2.32 i, s = 0.0416. (d) Flow profile shown in Figure 9 with 
M A v 9 = -0.115, best fit Z = 2.37 - 2.35 i, s = 0.0373. 
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(a) 1 Radial Basis Function 


Re{Z} 

(b) 2 Radial Basis Functions (c) 3 Radial Basis Functions 


Figure 8 . Contours of minimum least-squares fit error as a function of the wall impedance Z for 1, 2, or 3 radial basis 
functions (modes). Data from rotating rake measurements, m — — 8 , nominal wall impedance Z = 1.7 — 2.4i, 77 = 10.8, 
and M.Avg = —0.115 with mean flow profile shown in Figure 9. 


not known if this is due to the lower flow speed or that the lower flow speed profiles are flatter as shown in 
Figure 9 for the Mavq — 0.085 profile compared to the higher speed profile. More measurements at different 
frequencies and higher flow speeds are required. For these three in-duct, rotating rake measurements, the 
effect of the change in the flow profile, from plug flow to measured flow profile, and the effect of this flow on 
the boundary condition, equation (2), results in an increase in the measured wall resistance and an absolute 
increase in the measured wall reactance. Overall, the variability we see for the best fit impedances produced 
in this study is shown in Figure 10 to often be within the variability of the impedances measured by other 
techniques. 


IV. Discussion 

A numerical technique has been developed to compute the eigenvalues and eigen-functions for the duct 
solution with soft walls and a general mean flow. It uses the linear companion matrix method based on 
a discretization of the third-order governing equation and its boundary conditions. This method does not 
require an initial guess in order to solve the eigenvalue problem producing both the axial wave numbers, the 
eigenvalues, and their associated eigen-functions. The latter formed the basis functions that were used in the 
least-squares fit of the data. The results in this paper show that these basis functions can easily be used in a 
least-squares fit to determine the radial mode amplitudes as done in previous studies where closed form basis 
function solutions were used. With the numerically computed basis functions, the mean flow can be correctly 
accounted for whether it is a plug flow or a mean flow profile. It was shown in an example case using rotating 
rake measured data that when the mean flow profile is included, the mode amplitude increased 2.5 dB. This 
potential for changes in the mode amplitude has not yet been verified by measurements. It was mentioned 
above that once the mode amplitude is determined, the mode power can be calculated. Then assuming that 
the least attenuated modes are propagating out the duct, the in-duct mode power can be compared to far 
field radiated power. Thus, a check on the mode amplitude can be made with these measurements. This 
data could also provide a constraint on fitting the data to find the wall impedance. 

While radiating mode power is useful, an issue of importance is that the results for best fit impedance 
assume that the modes are propagating in the direction out of the duct which is the same assumption used 
when duct radiating power is computed. When in-duct, rotating rake measurements are made over acoustic 
treatment, modes may exist in the local vicinity that are propagating in both directions due to reflections in 
the duct. Currently, these are all lumped together in the measured data. The analysis cannot separate modes 
propagating in opposite directions. It only computes mode amplitudes that best fits the basis functions given 
and it can have free reign to adjust the impedance as necessary to further improve the fit. This could be a 
reason why the fit error versus impedance grid search has such a low variability in the size of the error leading 
to the situation in Figure 7 where multiple fit solutions were found for the same data. A way to separate 
the direction of propagating modes is to take measurements at separate axial locations simultaneously. 

Another issue to consider is that the expectation function for computing the residual sum of squares 
only used data from one circumferential m-order. Additional m-order equations can be included in equation 


NASA/TM— 2007-2 14966 


15 




(14) each with their own circumferential mode data from rotating rake measurements and with their radial 
basis functions individually computed at the same frequency, wall admittance, and mean flow profile. With 
multiple radial basis functions competing to fit multiple sets of measured data all with the same impedance 
boundary condition (neglecting any effects of mode amplitude on the impedance) , then perhaps the minimum 
error might be more definitive with errors away from minimum error being much larger than currently found 
in this study. Providing more data for the fit in that manner would be easier to do than, say, adding more 
microphones to the rotating rake. But, adding more microphones is still important for resolving higher 
orders of radial modes. 
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Figure 9. In-duct mean flow profiles. MAvg — 0.250 
profile derived from CFD calculations. MAvg = 0.115 
and 0.085 profiles show measurement points in the 
ANCF with lines for interpolated results. 


Figure 10. Wall impedance determined from best fit to 
Rotating Rake measured data. Open circle: plug flow 
calculations. Closed circle: flow profile calculations. 
All other data from Sutliff. 13 
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